Modélisation statistique

Préparation

Packages

library(dplyr)
library(ggplot2)
library(plotly)
library(glmnet)
library(gbm)

Chargement des données

load("../1_data_management_dplyr/fichiers_prepared.RData")

Jointure INSEE & FINESS

nb_ET_CODGEO=finess_et%>%
  group_by(CODGEO)%>%
  summarise(nb_ET=n())

sum(!nb_ET_CODGEO$CODGEO%in%insee$CODGEO)

## [1] 3191

nb_ET_CODGEO=merge(nb_ET_CODGEO,insee,by="CODGEO",all.y=T)
nb_ET_CODGEO <- nb_ET_CODGEO%>%mutate(nb_ET=ifelse(is.na(nb_ET),0,nb_ET))
table(nb_ET_CODGEO$nb_ET)%>%head

## 
##     0     1     2     3     4     5 
## 28153  2888  1400   947   583   425

Sélection des variables pour entraîner un modèle

data_model <-  nb_ET_CODGEO%>%
  select(-LIBGEO,-CODGEO)%>%
  mutate_if(is.character,as.numeric)

Echantillonage pour séparer apprentissage/test ie train / test

train=sample(x = nrow(data_model),size = round(.65*nrow(data_model)))

GLM

Préparation des données

Les valeurs manquantes sont très mal gérées par les GLM, on va donc imputer par la moyenne.
Il existe de nombreuses stratégies d’imputation donc certaines s’appuient sur des modèles de machine learning avancé : - Moyenne - Médiane - HotDeck - HotDeck stratifié - k plus proches voisins (kNN) - Régression/classification GLM - Régression/classification CART - Régression/classification GBM/RandomForest

Il existe plusieurs packages R qui implémente des stratégies d’imputation comme MICE ou simputation.
Lorsque les contraintes et le contexte sont très spécifiques et demande une maîtrise importante, il ne faut pas hésiter à implémenter son modèle d’imputation.
Vous pouvez vous référer à la présentation (Séminaire EIG #2) sur l’imputation des non-réponses partielles pour l’enquête OC : https://gitlab.com/DREES_code/OSAM/enquete_oc

data_model_imputed_for_glm=data_model%>%
  mutate_all(function(x)ifelse(is.na(x),mean(x,na.rm=T),x))

Entraînement du modèle

Un GLM dans R ce décrit comme une formule Y~X1+X2+X3 où Y est la variable à expliquer/prédire et X1, X2, X4 les variables explicatives.
Lorsqu’on a déjà sélectionné les variables pertinentes pour le modèle, il suffit d’écrit . pour utiliser toutes les variables disponibles sauf Y.
Si on a souhaité conserver la variable ID et qu’on souhaite garder toutes les variables sauf ID il suffit d’écrire Y~.-ID.

On peut ajuster le modèle avec un certain nombre de parmaètres facultatifs. - Choisir la fonction de lien (log, inverse, logit, identité…) et la loi conditionnelle. Si on a oublié quelle est la fonction de lien canonique associée à une loi de la famille exponentielle, pas de problème, elle est assignée par défaut. - On peut préciser un ou plusieurs offset (coeff fixé, pas forcément à 1). - Pondérer les observations avec weights (par exemple si on veut donner moins d’importance aux données plus anciennes). Des idées générales sur cette approche, valable pour la plupart des modèles prédictifs. ici. - On peut définir les contrasts pour préciser quelle catégorie doit être utilisée comme référence pour une variable explicative catégorielle.

model <- glm(nb_ET~.,
             data=data_model_imputed_for_glm[train,],
             family=poisson(link = "log"))

Coefficients du modèle

coeff_glm=summary(model)$coefficients
coeff_glm

##                      Estimate   Std. Error     z value      Pr(>|z|)
## (Intercept)     -1.857403e+01 1.112083e+00 -16.7020210  1.266995e-62
## NBMENFISC14     -2.202698e-04 3.130643e-06 -70.3592817  0.000000e+00
## NBPERSMENFISC14  1.080250e-04 1.536689e-06  70.2972213  0.000000e+00
## MED14            3.873753e-05 2.700597e-06  14.3440606  1.160611e-46
## PIMP14           5.959454e-02 3.305848e-03  18.0270052  1.195964e-72
## TP6014          -4.041158e-01 5.961887e-03 -67.7831956  0.000000e+00
## TP60AGE114      -1.648667e-03 1.774741e-03  -0.9289619  3.529088e-01
## TP60AGE214      -2.975953e-02 2.476067e-03 -12.0188701  2.828122e-33
## TP60AGE314      -2.221571e-02 2.829859e-03  -7.8504630  4.145040e-15
## TP60AGE414      -5.802431e-03 3.043607e-03  -1.9064328  5.659408e-02
## TP60AGE514       1.062246e-02 3.206787e-03   3.3124913  9.246901e-04
## TP60AGE614      -1.504853e-01 2.759352e-03 -54.5364658  0.000000e+00
## TP60TOL114      -9.450752e-03 4.351147e-03  -2.1720138  2.985462e-02
## TP60TOL214      -2.090915e-02 2.613153e-03  -8.0015040  1.229086e-15
## PRA14            1.879229e-01 8.423393e-03  22.3096460 2.978438e-110
## PTSAC14          4.259664e-02 5.906041e-03   7.2123841  5.498062e-13
## PPEN14           2.824090e-01 1.067305e-02  26.4600132 2.798566e-154
## PPAT14           2.905392e-01 1.029103e-02  28.2322845 2.348884e-175
## PPSOC14          7.946753e-01 1.110472e-01   7.1561937  8.294792e-13
## PPFAM14         -3.744029e-01 1.117036e-01  -3.3517524  8.030180e-04
## PPMINI14        -5.776451e-01 1.120897e-01  -5.1534198  2.557784e-07
## PPLOGT14         2.470051e+00 1.161699e-01  21.2624049 2.531083e-100
## D114            -7.712196e-04 2.328935e-05 -33.1146897 1.826512e-240
## D914            -1.493252e-04 5.239943e-06 -28.4974829 1.258660e-178
## RD14             1.611772e+00 5.182045e-02  31.1030042 2.193453e-212

Prediction

pred_glmtrain=data.frame(pred=predict(model,newdata = data_model_imputed_for_glm[train,],type="response"), obs=data_model[train,]$nb_ET)

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

perf_train=data.frame(pred=data_model[train,]$nb_ET, obs=data_model[train,]$nb_ET)


pred_glmtest=data.frame(pred=predict(model,newdata = data_model_imputed_for_glm[-train,],type="response"),obs=data_model[-train,]$nb_ET)

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

perf_test=data.frame(pred=data_model[-train,]$nb_ET, obs=data_model[-train,]$nb_ET)

Courbe de Lorenz ou Gain Curves

Il s’agit des courbes de mesure des inégalités de richesse.
Ici on les applique à la mesure des inégalités de Y (dans notre cas nb_ET) dans la population : - La courbe “perfection” décrit les vraies d’inégalités - La courbe “random” décrit un modèle incapable de discerner les inégalités - La courbe “glm” décrit la capacité de notre modèle à appréhender (car échantillon de test) les inégalités

L’intérêt d’une telle mesure est qu’on peut définir - ce qu’est un “mauvais” modèle (random) - ce qu’est un “bon” modèle (perfection) et comparer notre modèle à ces deux extrêmes.
Cette approche est inspirée de l’étude de la courbe de ROC pour les modèles binaires (Bernoulli).

L’inconvénient est que cette métrique mesure les inégalités en rang et pas en taille.

pred_glmtest%>%
    na.omit%>%
    arrange(-pred)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      data.frame(y100=quantile(.$y,0:100/100),
                 x100=quantile(.$x,0:100/100))
    }->Lorenz_glm_points
  
pred_glmtest%>%
    na.omit%>%
    arrange(-obs)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      data.frame(y100=quantile(.$y,0:100/100),
                 x100=quantile(.$x,0:100/100))
    }->Perfection
  
g <- ggplot()+
    geom_line(data = Lorenz_glm_points,
              aes(x=x100,y=y100,color="glm"))+
    geom_line(data = Perfection,
              aes(x=x100,y=y100,color="perfection")) +
    geom_line(data=data.frame(x100=c(0,1), yrandom=c(0,1)),
              aes(x=x100,y=yrandom,color="random"))
  
g 

Coeff de Gini

Par extension de l’AUC sous la ROC, on calcule l’AUC sous la courbe de Lorenz.
Cette métrique, translatée (X->2*X-1) pour se comparer à “random” est appelée indice de Gini.
Pour bien faire, on peut rapporter le Gini du modèle au Gini de la “perfection” pour que notre indice ait des valeurs POSSIBLES entre 0 et 1.

Gini=function(pred){
  pred%>%
    na.omit%>%
    arrange(-pred)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      # print(paste0("nombre_de_lignes:",nrow(.)," indice de Gini:"))
      2*mean(.$y)-1
      }
}
Gini(pred_glmtest)

## [1] 0.6163952

Gini(perf_test)

## [1] 0.9311837

Gini(pred_glmtrain)

## [1] 0.6401357

Gini(perf_train)

## [1] 0.9314016

Performance faible mais pas d’overfitting ie peu d’écart entre apprentissage et test.

GLMNET ou Elastic Net (mix L1 & L2)

Qu’est-ce que c’est ?

Très bonne explication ici. 2 paramètres à choisir : - alpha qui équilibre entre les pénalisations L1 et L2. - lambda qui équilibre la fonction de perte entre vraisemblance et pénalisation sur la taille des coefficients Lambda élevé signifie une pénalisation importante sur la taille des coefficients.
Lambda faible signifie qu’on est proche d’un GLM classique non pénalisé.

Le package R glmnet test automatiquement et de manière optimisé plusieurs lambda.
En revanche c’est à nous de tester plusieurs paramètres alpha. C’est ce qu’on appelle l’optimisation des hyperparamètres (hyperparameters tuning). La méthode qui consiste à tester les combinaisons des différents hyperparamètres s’appelle gridsearch.
La méthode gourmande/gloutonne (bruteforce) consiste à tester toutes les combinaisons.
Il existe des méthodes bayesiennes pour explorer intelligemment la grille de valeurs.
De nombreux packages R sont disponibles, mlrMBO est plutôt bien maintenu.
Mais ces considérations dépassent l’ambition de cette formation.

Préparation des données

Dans R il y a beaucoup de liberté quant au traitement des données avec les formats matrix, data.frame, data.table…
Pour glm, on avait un modèle de la forme cible ~ variables explicatives ou Y~X avec Y et X dans un même data.frame.
Pour glmnet, on doit fournir Y comme un vecteur et X comme une matrice de valeurs numériques.
Par conséquent les variables catégorielles devront être binarisées (dummy variables), ce procédé s’appelle le one-hot encoding. C’est très simple en R, il suffit d’utiliser la fonction model.matrix. Un exemple ici
Pas besoin d’y recourir dans notre cas. Si

target_train = data_model_imputed_for_glm[train,]$nb_ET
explanatory_train = data_model_imputed_for_glm[train,]%>%
  select(-nb_ET)%>%as.matrix

target_test = data_model_imputed_for_glm[-train,]$nb_ET
explanatory_test = data_model_imputed_for_glm[-train,]%>%
  select(-nb_ET)%>%as.matrix

Entraînement du modèle

thresh est le coefficient de convergence. On peut lire dans la documentation :
thresh Convergence threshold for coordinate descent. Each inner coordinate-descent
loop continues until the maximum change in the objective after any coefficient
update is less than thresh times the null deviance. Defaults value is 1E-7.

glmnet_model <- glmnet(x=explanatory_train,
       y=target_train,
       family="poisson",
       alpha=1,
       nlambda=100,#par défaut
       thresh = 1e-06,#par défaut
       maxit=1e7)
plot(glmnet_model$lambda,glmnet_model$dev.ratio)

# Distribution des coefficients en fonction du lambda.
plot(glmnet_model, xvar='lambda',label=T)

s=sample(glmnet_model$lambda,1)
smin=min(glmnet_model$lambda)
  pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=0)[,1],obs=target_test)
  print(paste("s=0",Gini(pred_glmnettest)))

## [1] "s=0 0.607603151181879"

  pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)
  print(paste("s=smin=",smin,"Gini=",Gini(pred_glmnettest)))

## [1] "s=smin= 0.000592615419036586 Gini= 0.607603151181879"

    pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=s)[,1],obs=target_test)
  print(paste("s=random",s,"Gini=",Gini(pred_glmnettest)))

## [1] "s=random 0.0185234740211396 Gini= 0.593625995260899"

[Avancé] cross-validation et choix du lambda optimal

Pour automatiser et rationnaliser le choix du lambda on fait de la validation croisée (cross-validation) ie on coupe le dataset en NFOLDS, disons en 10, et entraîne sur 90% vs validation sur 10% 10 fois.

NFOLDS=10
for (alpha in 0:10/10){
  print(alpha)
  
  tryCatch(rm(glmnet_cv,pred_glmnettest),warning = function(w) {
    print(paste("warning",w))
}, error = function(e) {
    print(paste("error",e))
})
  tryCatch({
glmnet_cv = cv.glmnet(x = explanatory_train, y = target_train, 
                              family = "poisson",#'gaussian', 
                              # Pénalisation L1 100%
                              alpha = alpha,
                              # On s'intéresse à la deviance - on suppose que la distribution conditionnelle suit une loi de Poisson
                              type.measure = "deviance",#'rmse'
                              # 5-fold cross-validation
                              nfolds = NFOLDS,
                              # valeurs élevée pour un entraînement plus rapide mais moins performant
                              thresh = 3e-5,
                              # On limite le nombre d'itérations
                              maxit = 1e6)
s1se=glmnet_cv$lambda.1se
smin=glmnet_cv$lambda.min
  pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=0)[,1],obs=target_test)
  print(paste("s=0",Gini(pred_glmnettest)))
  pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)
  print(paste("s=smin",Gini(pred_glmnettest)))
    pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=s1se)[,1],obs=target_test)
  print(paste("s=s1se",Gini(pred_glmnettest)))
  },warning = function(w) {
    print(paste("warning",w))
}, error = function(e) {
    print(paste("error",e))
})
}

## [1] 0
## [1] "warning simpleWarning in rm(glmnet_cv, pred_glmnettest): objet 'glmnet_cv' introuvable\n"
## [1] "s=0 0.55603861745672"
## [1] "s=smin 0.471635636426183"
## [1] "s=s1se -0.0810306213816748"
## [1] 0.1
## [1] "s=0 0.583863340143182"
## [1] "s=smin 0.581864126049507"
## [1] "s=s1se 0.574891153759821"
## [1] 0.2
## [1] "s=0 0.584880645708941"
## [1] "s=smin 0.584040178833034"
## [1] "s=s1se 0.575241615129584"
## [1] 0.3
## [1] "s=0 0.585850283170722"
## [1] "s=smin 0.584491320149618"
## [1] "s=s1se 0.57168251761927"
## [1] 0.4
## [1] "s=0 0.585168021660252"
## [1] "s=smin 0.583594171970502"
## [1] "s=s1se 0.570500609548991"
## [1] 0.5
## [1] "s=0 0.586404015902074"
## [1] "s=smin 0.585704976725855"
## [1] "s=s1se 0.573490703512848"
## [1] 0.6
## [1] "s=0 0.585778456623295"
## [1] "s=smin 0.585172256189612"
## [1] "s=s1se 0.569911026329369"
## [1] 0.7
## [1] "s=0 0.585120862815982"
## [1] "s=smin 0.584173549067466"
## [1] "s=s1se 0.567828649427626"
## [1] 0.8
## [1] "s=0 0.585781386610661"
## [1] "s=smin 0.585090279328811"
## [1] "s=s1se 0.55780532310059"
## [1] 0.9
## [1] "s=0 0.585819469470257"
## [1] "s=smin 0.585195521684529"
## [1] "s=s1se 0.569333421177145"
## [1] 1
## [1] "s=0 0.5858346984284"
## [1] "s=smin 0.584967554715142"
## [1] "s=s1se 0.570398813416507"

Gini(pred_glmtest)

## [1] 0.6163952

Comparaison pour différents lambda

plot(glmnet_cv)

Coefficients du modèle pour le “meilleur” lambda

Des idées pratiques en anglais ici
et théoriques en français ici

s1se=glmnet_cv$lambda.1se
smin=glmnet_cv$lambda.min
coeff_glmnet=coef(glmnet_model,s=s1se)
coeff_glmnet

## 27 x 1 sparse Matrix of class "dgCMatrix"
##                             1
## (Intercept)      1.300955e+01
## NBMENFISC14     -3.511324e-05
## NBPERSMENFISC14  1.739912e-05
## MED14            1.939810e-05
## PIMP14           7.064657e-02
## TP6014          -2.495468e-01
## TP60AGE114      -2.815586e-03
## TP60AGE214      -3.390306e-02
## TP60AGE314      -1.941505e-02
## TP60AGE414       .           
## TP60AGE514       9.446590e-03
## TP60AGE614      -1.557449e-01
## TP60TOL114      -1.705526e-02
## TP60TOL214      -5.546446e-02
## PRA14           -2.575534e-02
## PTSAC14          .           
## PBEN14          -5.715774e-02
## PPEN14           1.041573e-02
## PPAT14           .           
## PPSOC14          .           
## PPFAM14          1.434758e-01
## PPMINI14         .           
## PPLOGT14         2.264310e+00
## PIMPOT14        -1.993671e-01
## D114            -1.176157e-03
## D914            -2.283935e-05
## RD14             2.580683e-01

Comparaison des modèles GLM vs GLMNET

Comparaison des coefficients

# On prépare la table des coeffs de glmnet
coeff_glmnet=as.matrix(coeff_glmnet)
coeff_glmnet=data.frame(var=rownames(coeff_glmnet),coeff_glmnet=coeff_glmnet[,1])
coeff_glmnet$var=as.character(coeff_glmnet$var)

# On prépare la table des coeffs de glm
coeff_glm=data.frame(var=rownames(coeff_glm),coeff_glm=coeff_glm[,1])
coeff_glm$var=as.character(coeff_glm$var)

# On joint les deux tables et on les compare
comparaison_coeff=merge(coeff_glm,coeff_glmnet,by="var")
comparaison_coeff%>%
  mutate(ratio_glm_sur_glmnet=coeff_glm/coeff_glmnet)%>%
  as.tbl

## # A tibble: 25 x 4
##    var               coeff_glm coeff_glmnet ratio_glm_sur_glmnet
##    <chr>                 <dbl>        <dbl>                <dbl>
##  1 (Intercept)     -18.6         13.0                     -1.43 
##  2 D114             -0.000771    -0.00118                  0.656
##  3 D914             -0.000149    -0.0000228                6.54 
##  4 MED14             0.0000387    0.0000194                2.00 
##  5 NBMENFISC14      -0.000220    -0.0000351                6.27 
##  6 NBPERSMENFISC14   0.000108     0.0000174                6.21 
##  7 PIMP14            0.0596       0.0706                   0.844
##  8 PPAT14            0.291        0.                     Inf    
##  9 PPEN14            0.282        0.0104                  27.1  
## 10 PPFAM14          -0.374        0.143                   -2.61 
## # ... with 15 more rows

Corrélation entre les variables explicatives

Pour bien comprendre ce qui se passe, on jette un oeil à la matrice des corrélations. On se rend compte que lorsque deux variables sont très corrélées, en général le glmnet n’apprend que sur l’une des deux ie le coeff de l’autre passe à 0 ! C’est de la sélection de variables

cor(explanatory_train)%>%as.data.frame->cormat
rownames(cormat) <- colnames(cormat)
knitr::kable(cormat)

	NBMENFISC14	NBPERSMENFISC14	MED14	PIMP14	TP6014	TP60AGE114	TP60AGE214	TP60AGE314	TP60AGE414	TP60AGE514	TP60AGE614	TP60TOL114	TP60TOL214	PRA14	PTSAC14	PBEN14	PPEN14	PPAT14	PPSOC14	PPFAM14	PPMINI14	PPLOGT14	PIMPOT14	D114	D914	RD14
NBMENFISC14	1.0000000	0.9984167	0.0208183	0.0006392	0.0973090	-0.0898814	-0.0515755	0.0030850	0.0054769	-0.0258042	-0.1254197	-0.0541054	0.0089500	0.0423841	0.0384769	0.0089611	-0.0518556	0.0338667	0.0657956	-0.0179646	0.0795993	0.1032621	-0.1229715	-0.1317018	0.0986816	0.2568275
NBPERSMENFISC14	0.9984167	1.0000000	0.0221946	-0.0006170	0.1097336	-0.0857774	-0.0456387	0.0111994	0.0177155	-0.0107857	-0.1155511	-0.0419783	0.0165958	0.0544405	0.0524286	-0.0028750	-0.0676709	0.0229050	0.0795501	0.0018112	0.0887001	0.1134380	-0.1215908	-0.1393211	0.0959901	0.2609702
MED14	0.0208183	0.0221946	1.0000000	0.4126685	-0.2983518	-0.1527980	-0.2074498	-0.2348152	-0.1952208	-0.1325584	-0.0779240	-0.1663771	-0.2707437	0.2593144	0.2349316	0.0571125	-0.1534607	0.1623956	-0.3796705	-0.3243665	-0.3377208	-0.3673480	-0.3842982	0.3936242	0.4366865	0.1411240
PIMP14	0.0006392	-0.0006170	0.4126685	1.0000000	-0.6738356	-0.3535934	-0.4680754	-0.5191563	-0.4349506	-0.2984816	-0.1905955	-0.4106980	-0.6194433	0.5836097	0.5371525	0.0882601	-0.2731589	0.1879644	-0.8081711	-0.6106288	-0.7621252	-0.7947379	-0.8212543	0.8389193	0.6728852	0.0419388
TP6014	0.0973090	0.1097336	-0.2983518	-0.6738356	1.0000000	0.4358103	0.6234363	0.7300184	0.6314476	0.4888253	0.3046256	0.6110380	0.8183820	-0.3295723	-0.2935661	-0.0965909	0.0902070	-0.2267538	0.8498794	0.5583382	0.8633294	0.8253156	0.4404471	-0.7419392	-0.3851240	0.2382170
TP60AGE114	-0.0898814	-0.0857774	-0.1527980	-0.3535934	0.4358103	1.0000000	0.6802413	0.5688741	0.5992714	0.4387827	0.3177686	0.4159850	0.4789988	-0.2102842	-0.1978503	-0.0112030	0.0656244	-0.1342850	0.4646945	0.3594114	0.4518370	0.4279787	0.3155012	-0.2397163	-0.2660402	-0.0792504
TP60AGE214	-0.0515755	-0.0456387	-0.2074498	-0.4680754	0.6234363	0.6802413	1.0000000	0.8225507	0.7809385	0.5406213	0.2953815	0.5947279	0.6543600	-0.2743439	-0.2533915	-0.0372491	0.0885162	-0.1723208	0.6108186	0.4619510	0.6090881	0.5511701	0.3966656	-0.3641218	-0.3424759	-0.0496777
TP60AGE314	0.0030850	0.0111994	-0.2348152	-0.5191563	0.7300184	0.5688741	0.8225507	1.0000000	0.7979844	0.5580182	0.3329651	0.6622644	0.7551814	-0.2727380	-0.2466251	-0.0623077	0.0695433	-0.2030926	0.6754097	0.4926014	0.6806969	0.6170958	0.4231891	-0.4579585	-0.3616415	0.0218190
TP60AGE414	0.0054769	0.0177155	-0.1952208	-0.4349506	0.6314476	0.5992714	0.7809385	0.7979844	1.0000000	0.7334258	0.4556255	0.7313368	0.6189096	-0.1516406	-0.1318821	-0.0597123	-0.0403063	-0.2189225	0.5984320	0.4878360	0.6044549	0.4957400	0.3763708	-0.3621051	-0.3065579	0.0127077
TP60AGE514	-0.0258042	-0.0107857	-0.1325584	-0.2984816	0.4888253	0.4387827	0.5406213	0.5580182	0.7334258	1.0000000	0.6649183	0.7639708	0.4418234	-0.0113203	-0.0068105	-0.0189768	-0.1542809	-0.1892656	0.4522552	0.3958513	0.4890909	0.3018251	0.2722866	-0.2387701	-0.1972543	0.0328357
TP60AGE614	-0.1254197	-0.1155511	-0.0779240	-0.1905955	0.3046256	0.3177686	0.2953815	0.3329651	0.4556255	0.6649183	1.0000000	0.6638655	0.2766219	0.0147030	0.0081161	0.0281381	-0.1269813	-0.1414807	0.2937237	0.2289557	0.3894672	0.1198557	0.1943359	-0.1247814	-0.1158061	0.0189285
TP60TOL114	-0.0541054	-0.0419783	-0.1663771	-0.4106980	0.6110380	0.4159850	0.5947279	0.6622644	0.7313368	0.7639708	0.6638655	1.0000000	0.5915844	-0.1296592	-0.1296774	0.0298597	-0.0553549	-0.1688146	0.5304987	0.4201749	0.5981184	0.3609990	0.3280408	-0.3445819	-0.2347278	0.0740005
TP60TOL214	0.0089500	0.0165958	-0.2707437	-0.6194433	0.8183820	0.4789988	0.6543600	0.7551814	0.6189096	0.4418234	0.2766219	0.5915844	1.0000000	-0.3215933	-0.2891907	-0.0811822	0.1041602	-0.2323753	0.7304128	0.5483034	0.7254101	0.6680636	0.4948594	-0.5779092	-0.4286201	0.0373501
PRA14	0.0423841	0.0544405	0.2593144	0.5836097	-0.3295723	-0.2102842	-0.2743439	-0.2727380	-0.1516406	-0.0113203	0.0147030	-0.1296592	-0.3215933	1.0000000	0.9782652	-0.1256334	-0.8661113	-0.3454740	-0.3359036	0.0086203	-0.4181226	-0.4344144	-0.4395265	0.5105751	0.3955586	0.0314028
PTSAC14	0.0384769	0.0524286	0.2349316	0.5371525	-0.2935661	-0.1978503	-0.2533915	-0.2466251	-0.1318821	-0.0068105	0.0081161	-0.1296774	-0.2891907	0.9782652	1.0000000	-0.3286173	-0.8620158	-0.4062010	-0.2729723	0.0743442	-0.3718541	-0.3705829	-0.3466464	0.4714192	0.3258423	-0.0140542
PBEN14	0.0089611	-0.0028750	0.0571125	0.0882601	-0.0965909	-0.0112030	-0.0372491	-0.0623077	-0.0597123	-0.0189768	0.0281381	0.0298597	-0.0811822	-0.1256334	-0.3286173	1.0000000	0.1792827	0.3698694	-0.2239570	-0.3164273	-0.1253565	-0.2056437	-0.3434493	0.0700980	0.2427212	0.2102727
PPEN14	-0.0518556	-0.0676709	-0.1534607	-0.2731589	0.0902070	0.0656244	0.0885162	0.0695433	-0.0403063	-0.1542809	-0.1269813	-0.0553549	0.1041602	-0.8661113	-0.8620158	0.1792827	1.0000000	0.1717902	0.0708519	-0.2426255	0.1756882	0.1950789	0.1696866	-0.2796913	-0.2427807	-0.0638055
PPAT14	0.0338667	0.0229050	0.1623956	0.1879644	-0.2267538	-0.1342850	-0.1723208	-0.2030926	-0.2189225	-0.1892656	-0.1414807	-0.1688146	-0.2323753	-0.3454740	-0.4062010	0.3698694	0.1717902	1.0000000	-0.3987059	-0.5493964	-0.2759594	-0.3001348	-0.3562042	0.1745923	0.4472082	0.3380026
PPSOC14	0.0657956	0.0795501	-0.3796705	-0.8081711	0.8498794	0.4646945	0.6108186	0.6754097	0.5984320	0.4522552	0.2937237	0.5304987	0.7304128	-0.3359036	-0.2729723	-0.2239570	0.0708519	-0.3987059	1.0000000	0.7884788	0.9430116	0.9505255	0.6206068	-0.8066904	-0.6005253	0.0264911
PPFAM14	-0.0179646	0.0018112	-0.3243665	-0.6106288	0.5583382	0.3594114	0.4619510	0.4926014	0.4878360	0.3958513	0.2289557	0.4201749	0.5483034	0.0086203	0.0743442	-0.3164273	-0.2426255	-0.5493964	0.7884788	1.0000000	0.5703837	0.6417663	0.6343388	-0.5044596	-0.6610101	-0.3052755
PPMINI14	0.0795993	0.0887001	-0.3377208	-0.7621252	0.8633294	0.4518370	0.6090881	0.6806969	0.6044549	0.4890909	0.3894672	0.5981184	0.7254101	-0.4181226	-0.3718541	-0.1253565	0.1756882	-0.2759594	0.9430116	0.5703837	1.0000000	0.9015407	0.5171414	-0.7954651	-0.4681748	0.1729280
PPLOGT14	0.1032621	0.1134380	-0.3673480	-0.7947379	0.8253156	0.4279787	0.5511701	0.6170958	0.4957400	0.3018251	0.1198557	0.3609990	0.6680636	-0.4344144	-0.3705829	-0.2056437	0.1950789	-0.3001348	0.9505255	0.6417663	0.9015407	1.0000000	0.5527669	-0.8426726	-0.5366146	0.1199257
PIMPOT14	-0.1229715	-0.1215908	-0.3842982	-0.8212543	0.4404471	0.3155012	0.3966656	0.4231891	0.3763708	0.2722866	0.1943359	0.3280408	0.4948594	-0.4395265	-0.3466464	-0.3434493	0.1696866	-0.3562042	0.6206068	0.6343388	0.5171414	0.5527669	1.0000000	-0.5495386	-0.8325880	-0.4598725
D114	-0.1317018	-0.1393211	0.3936242	0.8389193	-0.7419392	-0.2397163	-0.3641218	-0.4579585	-0.3621051	-0.2387701	-0.1247814	-0.3445819	-0.5779092	0.5105751	0.4714192	0.0700980	-0.2796913	0.1745923	-0.8066904	-0.5044596	-0.7954651	-0.8426726	-0.5495386	1.0000000	0.5311730	-0.2764314
D914	0.0986816	0.0959901	0.4366865	0.6728852	-0.3851240	-0.2660402	-0.3424759	-0.3616415	-0.3065579	-0.1972543	-0.1158061	-0.2347278	-0.4286201	0.3955586	0.3258423	0.2427212	-0.2427807	0.4472082	-0.6005253	-0.6610101	-0.4681748	-0.5366146	-0.8325880	0.5311730	1.0000000	0.6519479
RD14	0.2568275	0.2609702	0.1411240	0.0419388	0.2382170	-0.0792504	-0.0496777	0.0218190	0.0127077	0.0328357	0.0189285	0.0740005	0.0373501	0.0314028	-0.0140542	0.2102727	-0.0638055	0.3380026	0.0264911	-0.3052755	0.1729280	0.1199257	-0.4598725	-0.2764314	0.6519479	1.0000000

Prediction du GLMNET

pred_glmnettrain=data.frame(pred=predict(glmnet_model,newx = explanatory_train,type="response")[,1],obs=target_train)

pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)

Courbe de Lorenz

pred_glmnettest%>%
    na.omit%>%
    arrange(-pred)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      data.frame(y100=quantile(.$y,0:100/100),
                 x100=quantile(.$x,0:100/100))
    }->Lorenz_glmnet_points


g <- g +
    geom_line(data = Lorenz_glmnet_points,
              aes(x=x100,y=y100,color="glmnet"))
g%>%ggplotly

Coeff de Gini

Gini(pred_glmtest)

## [1] 0.6163952

Gini(pred_glmnettest)

## [1] 0.607598

Gini(perf_test)

## [1] 0.9311837

Gini(pred_glmtrain)

## [1] 0.6401357

Gini(pred_glmnettrain)

## [1] -0.01661206

Gini(perf_train)

## [1] 0.9314016

GBM (gradient boosting method)

Préparation des données

L’avantage des modèles basés sur des arbres de décision, c’est qu’un arbre de décision sait bien gérer les NA.
Par exemple pour un arbre binaire, on peut imaginer séparer les NA dans un 3ème split, ou bien séparer NA vs le reste. La contrainte imposée par les variables continues est bien moins forte que dans un GLM. Avec un arbre si on coupe X à un seuil th (choisi par le modèle), on peut déciser de regrouper les NA avec le groupe 1 X<th ou avec le groupe 2 X>th.
L’imputation n’est pas nécessaire ici.

data_model%>%head

##   nb_ET NBMENFISC14 NBPERSMENFISC14   MED14  PIMP14  TP6014 TP60AGE114
## 1     0         304           799.5 21576.7      NA      NA         NA
## 2     0         100           235.5 21672.9      NA      NA         NA
## 3     0        6087         13660.5 19756.1 56.8215 15.7534    19.4181
## 4     0         603          1661.5 23204.8      NA      NA         NA
## 5     0          48           102.0 22157.5      NA      NA         NA
## 6     0        1030          2635.0 21679.3 61.7476  8.8425         NA
##   TP60AGE214 TP60AGE314 TP60AGE414 TP60AGE514 TP60AGE614 TP60TOL114
## 1         NA         NA         NA         NA         NA         NA
## 2         NA         NA         NA         NA         NA         NA
## 3    19.5204    19.1982    14.7159         NA         NA    5.40116
## 4         NA         NA         NA         NA         NA         NA
## 5         NA         NA         NA         NA         NA         NA
## 6         NA         NA         NA         NA         NA         NA
##   TP60TOL214 PRA14 PTSAC14 PBEN14 PPEN14 PPAT14 PPSOC14 PPFAM14 PPMINI14
## 1         NA    NA      NA     NA     NA     NA      NA      NA       NA
## 2         NA    NA      NA     NA     NA     NA      NA      NA       NA
## 3     24.796  71.8    67.9    3.9   27.3   10.1     6.5     2.8      1.8
## 4         NA    NA      NA     NA     NA     NA      NA      NA       NA
## 5         NA    NA      NA     NA     NA     NA      NA      NA       NA
## 6         NA  76.3    73.8    2.5   26.5    8.2     4.3     2.7      0.8
##   PPLOGT14 PIMPOT14    D114    D914    RD14
## 1       NA       NA      NA      NA      NA
## 2       NA       NA      NA      NA      NA
## 3      1.9    -15.7 10545.7 33376.5 3.16495
## 4       NA       NA      NA      NA      NA
## 5       NA       NA      NA      NA      NA
## 6      0.8    -15.3 12797.2 34818.9 2.72082

params=c(shrinkage=.01,
         nb_trees=1000,
         depth=10,
         nminobs=10,
         bag.frac=.2)

gbm_model <- gbm(nb_ET~.
                 ,data=data_model[train,]
                 ,verbose=T
                 ,train.fraction=.8
                 ,shrinkage = params[1]
                 ,n.trees = params[2]
                 ,interaction.depth = params[3]
                 ,n.minobsinnode = params[4]
                 ,bag.fraction = params[5]
                 )

## Distribution not specified, assuming gaussian ...
## Iter   TrainDeviance   ValidDeviance   StepSize   Improve
##      1      123.8052        156.1559     0.0100    1.1011
##      2      122.9204        155.2455     0.0100    0.9441
##      3      121.7714        153.9205     0.0100    1.1587
##      4      120.7843        152.6531     0.0100    0.9551
##      5      119.5231        151.0879     0.0100    1.1497
##      6      118.0686        149.2709     0.0100    1.3201
##      7      116.9232        147.9349     0.0100    1.1021
##      8      115.4956        145.9970     0.0100    1.2837
##      9      114.6430        145.0367     0.0100    0.8119
##     10      113.4236        143.6285     0.0100    1.1600
##     20      103.3187        131.8189     0.0100    0.7583
##     40       88.9986        116.3726     0.0100    0.1039
##     60       77.6587        103.1380     0.0100    0.4699
##     80       69.6424         94.2299     0.0100    0.2012
##    100       64.3673         88.7138     0.0100    0.2099
##    120       60.2482         84.8137     0.0100    0.2021
##    140       56.7311         81.3660     0.0100    0.0647
##    160       54.2119         79.3694     0.0100    0.0659
##    180       52.3049         77.7781     0.0100    0.0431
##    200       50.9097         77.1184     0.0100    0.0366
##    220       49.3800         75.8095     0.0100    0.0189
##    240       48.1694         74.8848     0.0100    0.0328
##    260       47.4105         74.7045     0.0100    0.0212
##    280       46.7316         74.6191     0.0100   -0.0054
##    300       46.1050         74.2809     0.0100    0.0045
##    320       45.0852         73.9441     0.0100   -0.0223
##    340       44.4359         73.7146     0.0100   -0.0258
##    360       43.5619         73.3486     0.0100   -0.0018
##    380       43.2005         73.3667     0.0100    0.0099
##    400       42.5862         73.5498     0.0100   -0.0096
##    420       41.9449         73.8726     0.0100   -0.0433
##    440       41.4313         73.6143     0.0100   -0.0003
##    460       40.7262         73.4490     0.0100   -0.0239
##    480       40.3588         73.5436     0.0100   -0.0230
##    500       40.0003         73.6608     0.0100    0.0022
##    520       39.4632         73.5346     0.0100    0.0038
##    540       39.0939         73.7524     0.0100    0.0080
##    560       38.7817         73.5695     0.0100   -0.0556
##    580       38.4370         73.3714     0.0100   -0.0259
##    600       38.1783         73.2017     0.0100   -0.0118
##    620       37.8702         73.3003     0.0100   -0.0118
##    640       37.4727         73.2190     0.0100   -0.0131
##    660       37.2582         73.3901     0.0100   -0.0208
##    680       36.8914         73.3819     0.0100   -0.0048
##    700       36.5775         73.6713     0.0100   -0.0281
##    720       36.3765         73.5775     0.0100   -0.0449
##    740       36.1241         73.5381     0.0100    0.0185
##    760       35.9153         74.0790     0.0100   -0.0128
##    780       35.6584         73.9118     0.0100   -0.0166
##    800       35.4054         74.1535     0.0100   -0.0450
##    820       35.0280         73.9798     0.0100   -0.0022
##    840       34.7990         74.2234     0.0100   -0.0183
##    860       34.4765         74.4036     0.0100   -0.0093
##    880       34.1561         74.8365     0.0100   -0.0102
##    900       33.7062         74.4357     0.0100   -0.0438
##    920       33.2471         74.2813     0.0100   -0.0056
##    940       32.9777         74.5498     0.0100   -0.0663
##    960       32.6899         74.3209     0.0100    0.0051
##    980       32.4597         74.1636     0.0100   -0.0072
##   1000       32.2198         73.8615     0.0100    0.0001

pred_gbmtrain=data.frame(pred=predict(gbm_model,
                                     newdata = data_model[train,]),
                        obs=target_train)

## Using 631 trees...

pred_gbmtest=data.frame(pred=predict(gbm_model,
                                     newdata = data_model[-train,]),
                        obs=target_test)

## Using 631 trees...

print(paste("gbm test:",Gini(pred_gbmtest)))

## [1] "gbm test: 0.846170486979527"

print(paste("gbm train:",Gini(pred_gbmtrain)))

## [1] "gbm train: 0.862547044056643"

print(paste("glm:",Gini(pred_glmnettest)))

## [1] "glm: 0.607597974870866"

Il y a plus d’overfitting (écart de 1 point entre train et test) dans le GBM que dans le GLM, mais il y a aussi un fit de bien meilleur qualité, on va comparer les courbes GLM et GBM sur le test.

Courbe de Lorenz

pred_gbmtest%>%
    na.omit%>%
    arrange(-pred)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      data.frame(y100=quantile(.$y,0:100/100),
                 x100=quantile(.$x,0:100/100))
    }->Lorenz_gbm_points


g <- g +
    geom_line(data = Lorenz_gbm_points,
              aes(x=x100,y=y100,color="gbm"))
g%>%ggplotly

Interprétation d’un GBM

Importance des variables

vars=summary(gbm_model)$var%>%as.character
summary(gbm_model)

##                             var    rel.inf
## NBMENFISC14         NBMENFISC14 31.1332462
## TP60AGE614           TP60AGE614 19.4038251
## NBPERSMENFISC14 NBPERSMENFISC14 16.2254886
## TP60AGE114           TP60AGE114  5.3784305
## RD14                       RD14  3.3676966
## TP60AGE514           TP60AGE514  3.0505128
## TP60TOL114           TP60TOL114  2.0291679
## PPAT14                   PPAT14  1.7369359
## PIMPOT14               PIMPOT14  1.6598815
## TP60AGE414           TP60AGE414  1.5248569
## PPLOGT14               PPLOGT14  1.3512653
## TP60AGE214           TP60AGE214  1.3377774
## PPEN14                   PPEN14  1.2507943
## TP60TOL214           TP60TOL214  1.1766011
## PBEN14                   PBEN14  1.1731137
## D914                       D914  1.1557584
## PIMP14                   PIMP14  0.9880151
## D114                       D114  0.9426558
## TP60AGE314           TP60AGE314  0.9134115
## PRA14                     PRA14  0.8606158
## PTSAC14                 PTSAC14  0.7240301
## PPMINI14               PPMINI14  0.6916395
## PPFAM14                 PPFAM14  0.6440630
## MED14                     MED14  0.5368247
## PPSOC14                 PPSOC14  0.3952722
## TP6014                   TP6014  0.3481201

Effet marginal moyen

Attention, contrairement à un GLM, ces courbes ne s’interprètent pas comme des effets toutes choses égales par ailleurs parce qu’elles s’appuient sur la distribution réelle des données.

plot(gbm_model,i.var=vars[1])

plot(gbm_model,i.var=vars[2])

plot(gbm_model,i.var=vars[3])

plot(gbm_model,i.var=vars[4])

A comparer avec les valeurs réellement prises par ces variables, les extrêmes sont estimés mais souvent les valeurs sont aberrantes et non significatives.
## Monotonie des variables

Si nécessaire on peut imposer la monotonie des variables, ceci réduira le surapprentissage et donnera des effets plus propres/plus interprétables.

dans gbm : var.monotone il convient de donner un vecteur nommé (avec les noms des variables explicatives) et des valeurs -1 si décroissante, 1 si croissante, 0 si absence de contrainte.

monotony=rep(0,ncol(data_model)-1)
names(monotony) <- setdiff(names(data_model),"nb_ET")
monotony["NBMENFISC14"] <- 1

params=c(shrinkage=.01,
         nb_trees=1000,
         depth=10,
         nminobs=10,
         bag.frac=.2)

gbm_model <- gbm(nb_ET~.
                 ,var.monotone = monotony
                 ,data=data_model[train,]
                 ,verbose=T
                 ,train.fraction=.8
                 ,shrinkage = params[1]
                 ,n.trees = params[2]
                 ,interaction.depth = params[3]
                 ,n.minobsinnode = params[4]
                 ,bag.fraction = params[5]
                 )

## Distribution not specified, assuming gaussian ...
## Iter   TrainDeviance   ValidDeviance   StepSize   Improve
##      1      123.9410        156.5011     0.0100    0.9179
##      2      122.8818        155.1943     0.0100    1.0843
##      3      121.9710        154.0022     0.0100    0.9235
##      4      120.6228        152.4411     0.0100    1.3434
##      5      119.5094        151.0565     0.0100    1.1113
##      6      118.2968        149.6567     0.0100    1.1069
##      7      117.0188        148.0968     0.0100    1.1690
##      8      115.9146        146.7008     0.0100    1.0543
##      9      114.6718        145.2054     0.0100    1.1690
##     10      113.7743        144.1456     0.0100    0.8839
##     20      103.6973        131.7029     0.0100    1.1680
##     40       89.2176        115.0668     0.0100    0.7386
##     60       78.3859        103.0026     0.0100    0.2495
##     80       71.1478         94.7057     0.0100    0.3714
##    100       65.4375         88.7413     0.0100    0.3448
##    120       60.7574         83.6699     0.0100    0.0602
##    140       57.5543         81.2544     0.0100    0.0992
##    160       54.8495         78.5961     0.0100    0.0180
##    180       52.4723         76.6235     0.0100    0.0541
##    200       50.7157         75.4813     0.0100    0.0454
##    220       49.4171         75.3343     0.0100   -0.0173
##    240       48.5624         74.7814     0.0100    0.0585
##    260       47.2398         74.1648     0.0100    0.0036
##    280       46.3777         73.8034     0.0100   -0.0079
##    300       45.3732         73.3304     0.0100   -0.0081
##    320       44.7860         73.5353     0.0100   -0.0355
##    340       44.1316         73.6588     0.0100   -0.0040
##    360       43.5558         73.4873     0.0100    0.0188
##    380       42.9367         73.7438     0.0100   -0.0111
##    400       42.3893         73.8192     0.0100   -0.0028
##    420       41.9995         73.7210     0.0100    0.0042
##    440       41.2510         73.8629     0.0100   -0.0126
##    460       40.8703         73.6792     0.0100   -0.0120
##    480       40.3314         73.4967     0.0100    0.0004
##    500       39.7779         73.1627     0.0100   -0.0690
##    520       39.3228         72.9770     0.0100   -0.0118
##    540       38.9252         72.7741     0.0100   -0.0202
##    560       38.5702         72.6034     0.0100    0.0128
##    580       38.2487         72.7149     0.0100   -0.0221
##    600       37.7219         72.4626     0.0100   -0.0137
##    620       37.4879         72.4712     0.0100    0.0066
##    640       37.1151         71.9383     0.0100   -0.0160
##    660       36.7128         72.4582     0.0100   -0.0360
##    680       36.4077         72.3271     0.0100   -0.0195
##    700       36.0229         72.0472     0.0100   -0.0025
##    720       35.7127         72.4194     0.0100    0.0017
##    740       35.2899         72.2755     0.0100   -0.0468
##    760       34.9286         72.4200     0.0100   -0.0209
##    780       34.6210         72.1668     0.0100   -0.0154
##    800       34.2715         71.5962     0.0100   -0.0263
##    820       33.9408         71.5788     0.0100   -0.0012
##    840       33.7161         71.7717     0.0100   -0.0260
##    860       33.3936         71.7323     0.0100   -0.0656
##    880       33.2113         71.2350     0.0100   -0.0174
##    900       33.0351         71.2736     0.0100   -0.0082
##    920       32.6330         71.2146     0.0100   -0.0292
##    940       32.4271         71.0239     0.0100   -0.0157
##    960       32.2729         71.2284     0.0100   -0.0067
##    980       31.8882         70.9893     0.0100   -0.0074
##   1000       31.6222         70.9821     0.0100   -0.0155

pred_gbmtrain=data.frame(pred=predict(gbm_model,
                                     newdata = data_model[train,]),
                        obs=target_train)

## Using 998 trees...

pred_gbmtest=data.frame(pred=predict(gbm_model,
                                     newdata = data_model[-train,]),
                        obs=target_test)

## Using 998 trees...

print(paste("gbm test:",Gini(pred_gbmtest)))

## [1] "gbm test: 0.844882059916177"

print(paste("gbm train:",Gini(pred_gbmtrain)))

## [1] "gbm train: 0.86362017133907"

print(paste("glm:",Gini(pred_glmnettest)))

## [1] "glm: 0.607597974870866"